Thermal Waves in Irradiated Protoplanetary Disks 

Sei-ichiro Watanabe 1 

Department of Earth and Planetary Sciences, Nagoya University, Chikusa, Nagoya 

464-8601, Japan 

seicoro@eps .nagoya-u. ac . jp 

and 
D. N. C. Lin 2 

UCO/Lick Observatory, University of California, Santa Cruz, CA 95064 

lin@ucolick . org 

ABSTRACT 

Protoplanetary disks are mainly heated by radiation from the central star. 
Since the incident stellar flux at any radius is sensitive to the disk structure 
near that location, an unstable feedback may be present. Previous investigations 
show that the disk would be stable to finite-amplitude temperature perturba- 
tions if the vertical height of optical surface is everywhere directly proportional 
to the gas scale height and if the intercepted fraction of stellar radiation is de- 
termined from the local grazing angle. We show that these assumptions may not 
be generally applicable. Instead, we calculate the quasi-static thermal evolution 
of irradiated disks by directly integrating the global optical depths to determine 
the optical surface and the total emitting area-filling factor of surface dust. We 
show that, in disks with modest mass accretion rates, thermal waves are sponta- 
neously and continually excited in the outer disk, propagate inward through the 
planet-forming domains, and dissipated at small radii where viscous dissipation 
is dominant. This state is quasi-periodic over several thermal timescales and its 
pattern does not depend on the details of the opacity law. The viscous dissi- 
pation resulting from higher mass accretion stabilizes this instability such that 
an approximately steady state is realized throughout the disk. In passive proto- 
stellar disks, especially transitional disks, these waves induce significant episodic 
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changes in SEDs, on the time scales of years to decades, because the midplane 
temperatures can vary by a factor of two between the exposed and shadowed 
regions. The transitory peaks and troughs in the potential vorticity distribution 
may also lead to baroclinic instability and excite turbulence in the planet-forming 
regions. 

Subject headings: accretion, accretion disks — circumstellar matter — instabili- 
ties — planetary systems: protoplanetary disks — stars: pre-main- sequence - 
solar system: formation 



Introduction 



It has become widely accepted that dusty protoplanetary disks are heated by radiation 
from the central star, and that this heating mainly determines the physical structure of the 
outer regions of these disks. Observed infrared spectral energy distributions (SEDs) of the 
T Tauri disks imply that their effective temperatures T decreases with disk radius r more 
slowly than T oc r~ 3 / 4 . This temperature distributions is usually explained by a model in 
which the thermal structure of the disk is assumed to be geometrically flared, i.e., the surface 
height z s where stellar r adiation is absorbed curves away from the midplane (or equivalently 
z s /r oc r 7 , with 7 > 0) (lAdams. Lada. fc Shulll987l ). 



The outer regions of these disks are irradiated by the central star (IKenyon fc Hartmann 



19871 ) and the flaring enables the disks to absorb more radiation from the central star. 
In a steady state, the flaring index 7 of a purely irradiated optically-thick disk can be 
obtained from the balance between intercepted stellar flux F s incident on the surface at a 
low angle 9 and emitted blackbody flux from the disk interior, under the assumptions that 
1) the surface height z s is everywhere proportional to the vertical gas scale height h, 2) the 
intercepted fraction of stellar radiation is the sine of the local grazing angle (F s oc sin#), and 
3) the central star is a point source. Under these circumstances, ther e is a self-consistent 
power-law solution with 7 = 2/7, which corresponds to T oc r~ 3 / 7 ( jKusaka et al.l Il970 
Chiang fc Goldreichlll997r ) . 



In addition to this particular power-law solution, there exists an one-parameter family 
of solutions for the disk structure, including the d iverging aspect-ra tio solutions and the 
asymptotically conical (in which z s oc r) solutions (iDullemondl l2000l ) . Small variations in 
the values of z s at inner radii where the i n tegra tion starts can cause large differences in the 
disk structure at large radii. IDullemondl (120001 ) speculated that such a sensitive nature of 
the steady solutions is suggestive to an intrinsic instability which must be analyzed with 
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time-dependent governing equations. 

More realistic steady solutions can be obtained numerically or semi-analytically to take 



the viscous dissipation associated with the mas s -accretion flow (e.g . . IChiang fc Go 



19971: IChiang et all l200ll ; iDullemond et all 12001 



. Chiang & Go 


dreich 


D'Alessio et al. 


2006; 



Garaud fc Linll2007l ). The most important novel feature of these series of second-generation 
models is the assumed presence of super heated surface dust layers above and below the 
disk midplane (IChiang &: Goldreichl 119971 ). Grains in these layers are directly exposed to 
the stellar flux. Grains much smaller than the peak wavelength of the self emission are 
superheated because of their low emissivity. The disk interior is heated by the superheated 
dust of the layers rather than directly by the central star. 

The two-layer disk model clearly explain the silicate an d water-ice emissio n bands in 
observed SEDs of Herbig Ae/Be stars and T Tauri stars (IChiang et aJj 1200 ll ). But, in 
spite of their triumph in the modeling of the observed SEDs, these model^j] are based on 
the assumption that z s is pro portional to the ga s scale height h, with a fixed constant of 
proportionality x — z B /h = 4 (IChiang et al.ll200ll ). While this assumption has been justified 
by estimates which suggest that changes of x is small throughout the disk, the amount of dust 
in the superheated layer is very sensitive function of x- This dependence arises because the 
dust spatial density at z s = \h is proportional to exp[— x 2 /2]. Thus, the fixed- x assumption 
may cause large discrepancy in determined values of z s . 

In previous analy ses, the magnitude F s is directly determined from the grazing angle 
9 (IChiang et al.ll200ll ). This approximation justified only in the case that the length of 
absorption layer (say, where the optical depth to the starlight changes from 0.1 to 1) along 
the starlight is smaller than the lengths of radial variations of surface density or temperature. 
This condition would not be satisfied if the disk surface contains fluctuations resulting from 
the growth of short-wavelength perturbations (see § [2]) . 

In principle, z s is determined by the condition that the visual optical depth, which can 
be obtained by a direct integration along the rays of starlight, is unity and the surface filling 
factor A s of the irradiated dust grains can be calculated through the vertical (in the direction 
normal to the disk plane) integration of geometrical opacity times mass density from z s to 
infinity. This global procedure yields z s at any given radius which depends not only on 
the local value of h but also the disk structure interior to that radius. Thus, x = z s/h 



1 Some author (e.g., Dullemond et al. 2001 ; Tanaka et al. 20051 : Garaud fc Lin 2007 ) determined \ se lf~ 
consistently, but these calculations are based on the grazing-angle approximation and the assumption that 
X changes slowly with radius. 
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varies both in space and time. Following this procedure, we can check for self consistency by 
recalculating z s based on the temperature distribution obtained from the steady, constant-x 
model. With this inductive approach, we demonstrate that there are substantial differences 
between the iterated values of z s and the initial, assumed values of \h. We also find that 
A s calculated from the deduced values of z s is substantially different from the values of sin 9 
extrapolated from the constant-x model. 

Since the irradiation heating can play such a major role in determining the vertical 
structure of protoplanetary disks, it is important to investigate the stability of such disks 
against the excitation of ripples on their surfaces. Under the assumption t hat the thermal 



timescale is much longer than the dynamical timescale in protostellar disks, iD'Alessio et al. 



(119991 ) investigated the thermal stability of the irradiation-dominated disks, using a simple 
cooling equation. They found the vertically isothermal disk to be stable against finite am- 
plitude perturbations. The initial temperature perturbations propagate inward and damp 
out quickly. However, in their analysis, they assumed that x is constant throughout the disk 
and F s is given by 9 as in the grazing angle approximation. The inferred stability in that 
study may depend on these assumptions. As we will see in § [2j changes in \ may lead to an 
instability. 



With a linear perturbation analysis, iDullemondl (120001 ) showed that the flaring disk so- 



lution may become unstable to infinitesimal hydrodynamic perturbations w hen the cooling 



time of the disk is much shorter than the dynamical time (the opposite limit of ID'Alessio et al. 



19991 ). The amplitude of inwardly propagating waves grow s exponentially with a r a te that 



is a decreasing function of the wavelength. Subsequently, iDullemond fc Dominiki (j2004al ) 
constructed a series of numerical models to examine the two-dimensional structure and evo- 
lution of protoplanetary disks around Herbig Ae/Be stars. In these simulations, they studied 
the radiative transfer process under the assumption that the disk always maintains a hydro- 
static equilibrium. (This assumption would not be appropriate for the limit that the time 
scale of cooling is shorter than that of dynamics.) They found two sets of asymptotically 
steady-state results which include the monotonically flaring solutions and the self-shadowed 
solutions. For the second set of solutions, the disk has a puffed-up inner rim. They stated 
that their iteration procedure, in which the hydrostatic equilibrium and the radiative trans- 
fer are treated separately in alternate steps, may not, under some circumstances, lead to 
a set of converged solutions. In their simulations, some wave-like disturbances appear to 
propagate over the disk from one iterative step to the next and these transitory features are 
never damped out completely. Although they limited their presentations on disks around 
Herbig Ae/Be stars, where the perturbation on the disk structure by these waves are rela- 
tively minor, they revealed that this problem appears to be more serious for disks around 
T Tauri stars. We speculate that this perturbation may be related to the above-mentioned 
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instabilities operating in the irradiation-dominated outer regions of protostellar disks. 

In this paper, we attempt to address two issues: 1) Are the steady-state solutions of irra- 
diated disks constructed by the previous one-dimensional models self-consistent and stable? 
2) Do these regions of disks tend to undergo quasi-periodic oscillations rather than attain 
an asymptotic steady state? In principle, these questions should be addressed with com- 
prehensive two- or three-dimensional numerical simulations. Such an approach is, however, 
fairly complicated, time consuming, and often plag ued with problems in the algor ithm which 



implements the radiative transfer processes (e.g.. iDullemond fc Dominikl l2004al ) . Prior to 
these detailed simulations, it is useful to identify the dominant effects which regulate the 
dynamics of irradiated disks with a set of one-dimensional time- dependent analyses on the 
thermal evolution of protostellar disks. 

Since thermal timescale is much shorter than viscous diffusion timescale, most previous 
studies adopted the assumption of thermal equilibrium during the course of disks' global 
evolution. Neverthele ss, there have been a few investigations on the thermal evolution of 



protoplanetary disks. IWatanabe et al.l (119901 ) investigated the cooling and quasi-static con- 
traction of the protoplanetary disks from an initial high-temperature state. They performed 
vertical one-dimensional numerical calculations and found that the cooling times are well 
estimated by a simple two-temperature (surface and interior temperatures) prescription. In 
this paper, we utilize this two-layer disk-temperature prescription to examine the stability 
and thermal evolution of irradiated disks. 

The simplest treatments for the thermal evolution of the irradiated disk are radial one- 
dimensional models in which the vertical structure of the disk at each radii is analyzed 
independently. In order to take into account of the irradiated surface and the disk interior, 
we evaluate the surface height directly from the location where the visual optical depth 
along the straight lines from the star is unity. We show that such disks evolve to quasi- 
periodic states in which thermal waves propagate inward through intermediate disk radii, 
where planets are formed. 

A simple discussion about the nature of thermal instability is given in § [2J The basic 
assumption of our model and its governing equations are presented in § O The results of our 
numerical calculations are presented in § H] for both simple and realistic opacities. Finally, 
we summarize our findings and discuss some possible evolutionary scenarios. 
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2. Simple discussion about thermal instability 

In this section, we discuss some potential causes for irradiation-dominated regions of 
disks to become thermally unstable. For illustration convenience, we adopt the following 
simplifying assumptions: 1) the star is a point source, 2) the disk's internal heat sources 
such as turbulent viscous heating as well as external heating other than the stellar radiation 
are negligible, 3) the optical depth of the disk is much larger than unity for both stellar 
radiation and its own emission, and 4) the transport of energy in the radial direction is much 
smaller than that in the vertical direction. Note that these assumptions are adopted only in 
this section for the purpose of pinpointing the physical process which leads to the thermal 
instability. All of these idealized assumptions will be relaxed in the numerical simulations 
to be presented below. 

Under these assumptions, we consider the heat balance in a geometrically thin disk 
which is irradiated by the central star. The energy equation reduces to 

BT 

C£^ = 2(F S -F m ), (1) 

where C is the specific heat per unit disk mass, E is the surface density of the disk, and 
T m is the temperature of the disk interior. Under assumption 2, the disk interior has an 
approximately isothermal structure and F m is the disk blackbody emission given by 

= aTl (2) 

and F s is the intercepted stellar flux given by 

where a is the Stefan-Boltzmann constant, L ± is the stellar luminosity, r is the cylindrical 
radial coordinate, and A s is the total emitting area-filling factor of superheated dust grains. 
The factor 1/2 in the right-hand side of equation ([3]) comes from the fact that surface 
irradiated dust re-radiate half of the absorbed stellar flux toward the disk interior (the rest 
toward infinity). Assuming a homogeneous mixing of gas and dust, we can obtain A s from 
(see Appendix [Cj 

A s = t g erfc ( J = r G erfc ( -= ) , (4) 



,V2hJ \V2 

where erfc(x) is the complimentary error function, tq is the geometrical optical depth of the 
disk midplane, and z s is the surface height where stellar radiation is absorbed. The ratio 
of the surface height z s to the gas scale height h is denoted by x = z s /h. In a hydrostatic 
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equilibrium, the gas scale height h is given by 



h n K \nm n GMj ' (5) 
where c m is the disk sound speed, Qk is the Keplerian angular velocity, k-Q is the Boltz- 
mann constant, \x is the molecular weight of disk gas, m u is the atomic mass unit, G is the 
gravitational constant, and M* is the stellar mass. 

Figure [1] displays the F m and F s (at 10 AU) as functions of T m . Assuming that the initial 
state is in a thermal equilibrium (corresponding to the point where three thick lines cross) 
with T m = T mieq , we impose a small positive temperature perturbation. We consider two 
extreme cases: 1) If z s is determined by the local disk structure, x — z s/h would retain a 
constant value during the increase of T m such that A s would also be constant (see eq. [I]) and 
F s would not change (dot- dashed line). In this case the system would be stabilized because 
F m > F s for T m > T m cq . 2) If z s is determined mostly by the attenuation by dust in the inner 
regions of the disk, z s would remain constant despite changes in the local disk temperature 
and h such that A s would increase rapidly and F s would increase much faster than F m (thick 
dash curve). In this case the system would be unstable because F m < F s for T m > T meq . 



In their stability analysis, iD'Alessio et al.l (119991 ) assumed a constant \- Based on 



the above analysis, this assumption naturally leads to stable solutions. In fact, most of 
the a nalysis on the stru cture of irradiated disks are based on the constant-x assumption 



Chiang et al.ll200ll ). The usual procedure to determine z s is based on a geometrical 



(e.g., 
consideration, i.e, 

^■ta^f^-lY (6 ) 
r \amr j 

where 9 is the grazing angle (i.e., the angle between the starlight and the disk surface). 
We refer this procedure to be the grazing-angle approximation. Most previous steady-state 
disk models are constructed with equation (jSJ) under the assumption that x is constant 
throughout the disk. 

However, the magnitude of \ is generally determined by the radial structure of the disk 
as well as its local properties. The surface height z s is determined by the optical depth 
integrated through a ray of the stellar radiation. We derive a ray integral and check the 
validity of equation (JHD in § [3] and Appendix ICl In principle, equations @ and (jSJ) must 



be resolved simultaneously (ITanaka et al.l 120051 ). However, this set of equations is fairly 



unstable to solve numerically because they do not contain contributions which may reduce 
any steep temperature gradients in the radial direction. 

Steep temperature gradient, if present, would invalidate the constant x an d the grazing- 
angle approximations. Physically, the radial transport of heat suppresses the radial temper- 
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ature gradient, but such a process through the opaque regions of the disk must be analyzed 
with multi-dimensional numerical simulations. One of the most efficient process of the radial 
heat transport is the radiative transfer from the superheated dust grains at the surface of 
any radial location to the disk midplane at adjacent radial regions. Using a simple one- 
dimensional model, we can take this oblique radiative transfer of heat into account. 



3. Basic Equations 



Following the approaches of IChiang fc Goldreichl (119971 ) and iGaraud fc Linl (120071 ) , we 
construct numerical models to study the thermal evolution of a protostellar accretion disk. 
The surface of the disk is illuminated by the central star. Exposed to the stellar radiation, 
sub-mm dust grains in the surface layers of the disk are superheated. We consider the case 
that dust mass of the disk is so large that the disk midplane (except for an innermost region 
where silicates are evaporated) is optically thick to the stellar radiation. In contrast to the 
previous section, the heat sources for the disk interior in these numerical models include 
both irradiation from the superheated grains on the disk surface and the viscous dissipation 
associated with the accretion flow. We adopt a cylindrical coordinate system (r, <f), z) in 
which the z = plane represents the disk midplane and the origin is at the location of the 
central star. Since the star-and-disk system is symmetric with respect to the midplane, we 
describe our results for the upper half of the disk only. 

In or der to simplify the p r oblem we adopt the two-layer axisymmetric disk model pro- 
posed by IChiang fc Goldreichl (119971 ). In this model, the disk consists of a surface super- 
heated layer where the dust temperature is T s (r) and a disk interior where the dust and gas 
temperature is assumed to be uniform at T m (r). This model is simple to use and includes 
all the essential ingredient to analyze the onset, evolution, and stabilization of thermal in- 
stability in protostellar disks. However, such a simplification would be invalid if the disk 
optical depth r m (T m ) to its intrinsic radiation is much larger than unity and viscous heating 
rate is larger than surface heating rate. However, the dust optical depth may be self-limited 
by their rapi d growth through co hesive collisions, so that r m (T m ) < 10 throughout the disk 
(see Fig. 6 in lTanaka et al.ll2005l ). Thus, the two-layer model is valid even in the inner disk 
where dust surface density is larger. 

The two-layer model is invalid at the inner edge of the disk, where disk is irradiated 
not only f rom the top but from the radi al direction. The disk may have a puffed-up inner 



rim (e.g., iDullemond fc Dominikl l2004al ) . but a set of two-dimensional radiative transfer 
calculations is needed to determine the structure of the innermost region. In this work we 
simply assume that disk within 0.1 AU is optically thick in the radial direction and has no 
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puffed-up rim that might cast shadows over the outer regions of the disk. We confine our 
calculations only in the regions r > 0.1 AU, where the two-layer model is valid. 

The thermal timescale of the disk interior at radius r is given by (see eq. |12j ) 

(7a + 1) cjE 
2( 7a - 1) oTi 



where 7 a is the adiabatic exponents and c m is the sound speed of the disk interior. For 
evaluation, we assume here power-law distributions for the total (gas + dust) surface density 
E(r) oc r~ p and the midplane temperature T m (r) oc r~ q . The normalization factors, So 
and T m0 , refer to their corresponding values at r = 1AU. The nominal value of surface 
density is given by that in the minimum mass solar-nebula (MMSN) model in which Sh = 



1.7 x 10 3 gcm" 2 flHavashilll981l ). 



We assume that the thermal timescale £ t h to be much longer than the dynamical time 
(fi^ 1 ), but much shorter than the viscous evolution time {r 2 /u, where v is the turbulent 
viscosity). In this case we can regard that the whole region of the disk is always in a 
hydrostatic equilibrium in the vertical direction and has time-independent surface densities. 

The temperature T s of superheated dust grains is given by 



4nr 2 



where e s is the averaged emissivity of the dust grains at T s . Along a ray from the surface of 
the star, the superheated layer extends outward until the position where visual optical depth 
has reached unity. We take into account the attenuation of the stellar photons by defining 
the height z s of the bottom of the superheated layer with the following equation 

r s (T ir ;r,z s (r)) = l. (9) 

Here T s (T*;r,z) is the optical depth between the central star and the point (r,z) to the 
blackbody radiation peaked at the stellar effective temperature T*, given by the following 
integration 

T B (T*;r,z) = f R B (n)p d (r',(r') (l + ( 2 f 2 dr> (10) 

along a straight path from the star to the point (r, z). Here R* is the stellar radius, £ = z/r 
is the aspect ratio, Pd(r', z') is the spatial mass density of dust at (r', z'), and /c s (T ra d) is the 
Planck mean opacity of the grains interacting with the blackbody radiation peaked at T rac j. 



-10- 



Note that we define here the grain opacity per unit dust mass, not per unit total (gas + dust) 
mass as in the usual definition, because it is convenient for the consideration of the case that 
the dust-to-gas ratio may change vertically. The emissivity in equation ([S]) can be given by 



The energy equation includes heating from both stellar irradiation an d viscous dissipa- 



tion as well as radiative losses from the disk surface such that (see, e.g., IWatanabe et al. 



19901) 



(7a + 1) k B EdT n 



2 [F a - F m ] + — Mfi|, 

47T 



2( 7a -l)/im u dt " rH ^ mJ ' 4*™-*' (12) 

where M is the mass accretion rate, which we assume to be constant throughout the disk. 
Note that the steady state assumption is compati ble with a power-law surfa c e density distri 



bution for some effective viscosity prescriptions (IChiang fc Goldreichl 119971 ; iD'Alessio et al. 
20061 : iGaraud fc LirJl2007h . 



Further, F s and F m are, respectively, the thermal radiation fluxes downward from the 
superheated dust grains high up in the disk atmosphere and upward from dust grains in the 
disk interior, 

U I A a AIL 



F s {r)= [1-e-^] 



in 



" 2 ^ 3nr 3 



(13) 



F m (r) = [l - e-*-PW] aTl (14) 

where r m (T s ) and r m (T m ) are the optical depths of the disk interior (from z = to z = z s ) to 
the radiation from the superheated dust grains and to its own emission, respectively, and A s 
is the total emitting-area filling-factor of superheated dust grains. We consider the effect of 
finite radius of the central star in F s , which is important in the inner part of the disk. We also 
consider the effects of oblique radiative transfer: the angular brackets in the right-hand side of 
equation ffl3l) represent radial averaging of radiation emitting from superheated dust within 
the adjacent regions (see §HJ). The factors 2 in the exponential fu nctions in equations ([TBI 
and (|14p also denote oblique radiative transfer in the disk interior (ITanaka et all 120051 ) . 



Once the dust density distribution p d (r,z) of the disk is specified, A s and r m can be 
determined from the following integration: 



A B (r) = 1 — exp 



H s {T*)pd(r, z ) dz 



z s (r) 



z a (r) 



K m (T iad )p d (r, z') dz 



(15) 
(16) 
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where « m (T rac i) is the Planck mean opacity of mi dp lane grains interacting the blackbody 
radiation peaked at temperature T rad . Further details of dust opacity are shown in §@] and 
Appendix [A] 

We assume that the total (gas + dust) surface density E of the disk is a simple power-law 
distribution in the radial direction: 



which is kept constant during the thermal evolution considered in this work. 

Taking the effects of dust settling, we can obtain dust density distribution p d of the 
disk in the two-temperature model adopted here (see Appendix [B]) . However, computation- 
ally intense iterative calculation is needed to determine, self-consistently, the magnitudes 
of pd and z s simultaneously. Instead of equations (1B4j) and (1B6I) . we adopt, in most of the 
calculations, the following simple density distribution for the dust: 



where Ed is the surface density of dust. We put Ed = /dE, where fd is the dust fraction in 
the surface density. We dub / d the dust-to-gas ratio. In the disk interior, dust sedimentation 
is not so important unless the radii of dust are not so large that we can formally set the dust 
density distribution to be equation ffTgj) . In the surface layer, dust settling reduces the dust 
density while high surface temperature T s raise it, so that equation (fl8l) also gives a good 
estimate. In some cases we compared results with more realistic dust distribution given in 
Appendix [B] and found that they are very similar unless dust sizes are not so large. We 
discuss the difference of the results between the two distributions in § 

We found that if the gradient dlnz s /d\nr changes rapidly in r direction, the approxi- 
mation used to derive equation is no longer justified (see Appendix [C]) . For this reason 
we use equations Q) and ffT5]) instead of equation (Q. 



We adopt the following values as fixed parameters for all models: the mass, radius, and 
luminosity of the central star are set to be M* = 1M , R* = 2.085-R Q , and L+ = 1L , 
respectively, so that its effective temperature is T* = 4000 K. For the disk gas, we specify 
p = 2.34 and 7 a = 1.4. 




(17) 




(18) 



4. 



Numerical Results 
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We adopted the phenomenological MMSN model (IHayashil Il98lh for the standard gas 
and dust surface density distribution. In this model, S = Xho = 1-7 x 10 3 gcm~ 2 with 
r = 1 AU and p = 1.5 in equation (IT7j) . For comparison, we also calculate a relatively flat 
E distribution with p = 1.0 and S = 3.54 x 10 2 gcm -2 . These surface densities are kept 
constant with time. The dust-to-gas ratio fd and solid material density p mat are set to be 
0.014 and 1.4 gem -3 , respectively. In § |4.1l we neglect any changes of dust surface density due 
to sublimation and use the value of /d throughout the disk. In § 14.21 we consider the effect of 
ice sublimation. We also vary the mass accretion rate from M = to lO~ 7 M yr _1 . These 
steady-state accretion rates are consistent with our specified surface density and temperature 
distributions provided the magnitude of a is a function of the radius, where a is the no n- 



dimensional turbulent viscosity in the so-called a-prescription (IShakura k, Sunyaevlll973l ) 



The normalization time unit is the thermal timescale given by equation (j7J) with S = S 
and T m = T m0 . We denote the time unit as £ t h,o an d the non-dimensional time as t = t/t t h,o- 
For the standard model, we set the magnitude of So = Sho and To = 124 K, so that 
£th,o — 53 yr. Note that the local thermal timescale t t t is nearly constant with r in the 
standard disk model with p — 1.5. The non-dimensional time step used in the numerical 
integration is set to be St = 0.005. In order to verify numerical convergence, we also perform 
several calculations with time step half of the standard value and confirmed that the results 
have no significant changes. 

The spatial grid consists of 90 points (the standard case) or 180 points (the high- 
resolution case), logarithmically distributed, between r = 0.1 AU and 100 AU. Numerical 
oscillation would be induced if there is no radial exchange of energy at all. At any radius, 
the disk interior is exposed not only to the superheated surface grains directly overhead, but 
also obliquely to those at adjacent radial locations. Such a radial exchange of energy tends 
to suppresses instabilities for short-wavelength perturbations. In our numerical scheme, we 
assume that the isotropic radiation comes from all radius r' within \r' — r\ < z s (r) contributes 
to the heating at radius r as expressed in equation ffl3l) . This implementation stabilizes the 
short wave-length oscillation and the results are essentially independent of the numerical 
resolution. 



4.1. Constant Opacity 



In this subsection, we first illustrate the dominant features using a simple opacity model. 
We adopt the emissivity and opacity (per unit dust mass) of the grains interacting with 
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blackbody radiation peaked at temperature Tj to be 

e s (Ti) = (Jpj and K B (Ti) = R s0 H±\ , (19) 

where we choose the value of R s q = 10 2 cm 2 g _1 , which approximately corresponds to the 
commonly defined opacity per unit gas mass with the value of lcm 2 g -1 . Most of the 
calculation shown here is for (3 = but we also calculate some models with /3 = 1 for 
comparison purposes. 

For initial conditions, we adopt the steady state solution obtained from the time inte- 
gration with fixing x{ r )- This set of initial conditions does not correspond to the asymptotic 
steady-state solutions because the initial estimate of x is n °t self-consistently compatible 
with the actual aspect ratio ( s of the surface. Nevertheless, the numerical calculations relax 
to nearby steady solutions if they exist. In order to verify that our results are independent 
of the adopted initial conditions, we calculate the evolution of the disk with several different 
initial guesses for x an d found that the system reaches to the same asymptotic state. 

We first show the results of the calculations with no mass accretion (M = 0). Figure [2] 
shows the initial evolution of the midplane temperature T m as well as the surface temperature 
T s , which is kept constant with time. The elapsed time is t — 1.6, i.e., t = 1.6t t h,o — 85 yr 
and each curve corresponds to a time step At = 0.2 (At = 10.6 yr). At first, the initial state 
is almost stable in the innermost region and the outermost region. But, in the intermediate 
regions (0.5-20 AU), the disk becomes unstable. At a typical instance of time, four local 
temperature peaks (high T m and z s ) coexist and are amplified. These peaks move inward 
(toward the star) as they grow. Hence, we refer these propagating perturbations as waves. 
During the amplification of the waves, they cast shadows over the outer regions of the peak. 
The temperature T m decreases in the shadowed regions. Each fully-grown wave has a sharp 
slope on the inner 'exposed side' of the peak and a gentler decline on the outer 'shadow side'. 
The waves propagate inward with velocities about a few tenth of r/t t h- The outermost wave 
(~ 16 AU) begins to grow just outside the shadowed region of the inner adjacent wave when 
the shadowed region is developed. This tendency shows that the outermost wave may be 
induced by the wave ahead of it. 

Additional time integration shows that as these waves propagate inward they begin 
to decay when they reach inside 1 AU. The waves are completely damped out at around 
0.25 AU. In contrast, new waves are formed continually at the outermost region (> 20 AU) 
of the computational domain. The growth region of waves gradually retreats and the maxi- 
mum amplitude of each wave during its propagation cycle gradually increases with time. The 
system reaches a quasi-periodic state when t ~ 8. Figure [3] shows the evolution of T m at this 
stage. Waves are continuously formed and amplified at the outer disk (> 30 AU), then prop- 
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agate inward with nearly constant amplitude through the intermediate disk regions, begin 
to decay at about 1 AU, and are damped out completely at around 0.25 AU. Temperature 
in the innermost disk region (r < 0.25 AU) attains steady values. The propagation speed 
of the waves is approximately given by r/t t h- m the intermediate disk radii (1-20 AU), the 
peak temperature of each wave is 2-3 times higher than the bottom temperature in the inner 
adjacent shadowed region. 

The radial profile of each wave is somewhat skewed. The half width of individual 
waves is about 0.1-0.2r on the inner side and 0.2-0.4r on the outer side. The wavelength 
is approximately twice as large as z 8 (r). Due to the steep radial temperature gradient, the 
magnitude of z s (r) is affected by the variation in the thickness h at the disk regions interior 
to r on this length scale. The ratio of the radii between two adjacent wave peaks are about 
20 if both waves are within r > 1 AU. Near the inner boundary of the propagating-wave 
zone, the finite size of the star (R*) becomes comparable to the surface height z s . This 
time-independent contribution in equation ({TBI essentially stabilizes the innermost region of 
the disk. 

The change of other variables at the same epoch as Figure [3] are shown in Figures 0HHJ 
Figure H] shows the time evolution of ( s = z s /r. Outside 0.25 AU, ( s changes stepwise. The 
two or three steep jumps of ( s correspond to the leading inner side of the thermal waves. The 
magnitude of increase in z s at the leading edge of each step reaches a maximum of 1.5 around 
10 AU. The flat portions of £ s correspond to the shadowed regions, where ( s is determined 
by the stellar ray which passes above the peak of each wave. 

Figure [5] shows the time evolution of x — z s/h = Cs/Ch- The local minima of x are 
located just ahead of the peaks in T m . The decrease of x is essential for the temperature 
raise. In the shadowed regions x increases because ( s is kept almost constant (see Fig. @J, 
whereas £h = h/r decreases. Changes in the value of x produces the variations in the 
surface filling factor A s . The time evolution of A s is shown Fig. [61 Sharp peaks of A s , which 
correspond to the minima of x ( see ec L- |Clj ). propagate inward. Note that the variation in 
the amplitude of A s is very large, ranging from a few tenth at the peaks to below 10~ 3 in 
the shadowed region. Such large changes of A s induce rapid heating in front of the waves 
and rapid cooling in the shadowed regions. The unperturbed z s in the intermediate regions 
is comparable to the half width of the propagating waves. Modest variations in x can l ea d 
to nonlinear dissipation, such that the wave amplitudes are also limited in these region. 

We also plot the evolution of the distribution of the logarithmic pressure gradient 
dlnP/dlnr (Fig. [7]), where P{r) is the pressure in the midplane of the disk. This plot 
indicates that the gas pressure gradient is nearly reversed just in front of the peak of the 
waves. This inversion occurs due to the steep positive temperature gradient in the exposed, 
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leading, inner face of the waves. Consequently, the velocity of the gas departs significantly 
from its unperturbed sub-Keplerian values. In a follow-up paper, we will consider the asso- 
ciated gas drag on grains of various sizes. Another interesting quantity is the distribution of 
the potential vorticity (or vortensity): 

^ = —Ar A n\ (20) 
S EH dr 

where Q ep is the epicyclic frequency and Q is the angular velocity of gas. Note that in a 
quasi-Keplerian disk Q ep ~ f2x oc r~ L5 , so that the potential vorticity is almost constant for 
MMSN with p = 1.5. The thermal waves disturb the potential vorticity through the change 
of pressure gradient. Figure [8] displays the evolution of the distribution of the potential 
vorticity. There are peaks and troughs around the waves. Local ex trema of this quantity 



can l ead baroclinic instabilities which may excite turbulence (e.g. iKlahr &: Bodenheimer 



2003 1) in the de ad zone where the magneto-rotational instability may have limited influences 



(lGammidll996l ). Further investigation of this possibility will also be considered elsewhere. 



Next, we examine the dependence of wave excitation and propagation on the mass 
accretion rate (M). Other than the value of M in the energy equation, we adopt the 
same model parameters as for the no-accretion case shown in Figures [2H5J Thus, in the 
present context, the primary physical effects associated with the accretion flow is the viscous 
dissipation. This internal energy source (viscous dissipation) in equation ( jl~2l) has a greater 
fractiona l contribution to the energy budget in the inner regions than the outer regions of 



the disk flGaraud & Linl 120071 ). 



Although quasi-periodic oscillations are excited in all cases, the radial extent where 
they propagate to depends on the magnitude of M. Figure M shows the time evolution of 
T m in a quasi-periodic state for the case of M = lO _8 M0yr _1 . The result is similar to that 
obtained by neglecting the viscous dissipation. In this case, waves are excited in the outer 
region, propagate inward and are damped out in the innermost region. In comparison with 
the no-accretion model, the innermost region in this case is hotter (T m (0.1AU) ~ 600 K) 
with larger thickness. The time-independent contribution of viscous heating to the energy 
equation provides a stabilizing effect to a slightly larger radial extent (r < 0.5 AU) compared 
with the no-accretion case (T m (0.1 AU) ~ 500 K and r < 0.25 AU). 

The quasi-periodic state, however, is drastically changed in the M = 1 x 10~ 7 M Q yr~ 1 
model (Fig. [TU]) . The disk becomes stable interior to about 6AU. In the outer region, there 
are two high-temperature peaks, both oscillate quasi-periodically. The positions of the two 
peaks do not coherently propagate inward as in the case of M — 1 x lO~ 8 M yr _1 , but 
they fluctuate to and fro around 12 AU and 40 AU, respectively. The ranges of temporal 
temperature changes are less than 2 K over the disk, so that the disk is regarded to be in 
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an approximately steady state. This result shows that the disk becomes stabilized as M 
increases. 

Finally, we show the dependence on the disk surface density S. Note that E affects the 
evolution not only through the optical depth of the surface layer but through the thermal 
timescale. We perform several calculations for So = 3.54 x 10 2 g cm -2 and p — 1.0. Note that 
the time unit for this case is tth,o = 11-0 yr. We calculate the disk evolution for this surface 
density distribution with several values of M and find that the results are quite similar to 
the case with the standard surface density distribution. The quasi-periodic wave solutions 
are obtained for small M, whereas the disk becomes nearly steady for large M. 

Figure [TT] shows the temporal variation of T m in the quasi-periodic state for this pre- 
scribed S(r) distribution with M = lO _8 M yr~ 1 . Comparing to the case with p = 1.5 
(Fig. ED, the wave propagating region is relatively narrow, confined to the regions between 
1 AU and 20 AU, and the maximum amplitude of waves is also smaller. Outermost region 
where r > 20 AU, the amplitudes of waves are limited and do not exceed about 2 K. The 
amplitudes grow as waves propagate from 20 AU to 8 AU, then they attain a nearly constant 
value from 8 AU to 1.5 AU. The waves are finally damped at 1.5-1 AU. In an analogous 
inviscid model (i.e., with an identically prescribed S(r) distribution but without viscous 
dissipation), the result is quite similar to that in Figure [TT1 except that the wave propagat- 
ing region extends slightly closer to the star. The results for the M = 10 _7 M Q yr -1 model 
show that the disk attains an approximately steady state with only very small fluctuations 
(amplitude is less than a few tenth K) remaining in the outermost region. 

All the above results are for (3 = 0. We also calculate cases with (3=1. Other than a 
modification in the distribution of T s and T m , the time dependent nature of wave excitation 
and propagation is essentially independent of the value of f3. In all cases, we find that 
the quasi-periodic nature of the inwardly-propagating thermal waves is realized for the disk 
with M < 10 _8 M Q yr _1 . The basic features of this state do not strongly depend on other 
parameters such as p and (3. In order to verify the universality of these results, we perform 
further calculations with more realistic opacities. 



4.2. Realistic Opacities 



In this subsection, we p resent the results ba sed on models with more realistic grain- 
opacity prescription given by iTanaka et al.l (120051 ) . In this prescription, grains are assumed 
to be consisted of a uniform mixture of H^O-ice, organics, olivine, pyr oxene, metallic iron, 
and troilite, the abundances of which are given by iPollack et al.l (I1994J ) . For relatively high 
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temperatures (T > 160 K), we use the dust opacity of grains without ice or organics. For 
the low-temperature (T < 160 K) state, the opacity includes the contribution from ice and 
organics grains. Ba sed on the single-sized monochromatic opacity table made by H. Tanaka 



( ITanaka et al.ll2005l ). we calculate the Planck mean of the size-averaged opacities K s (T rad ) and 
^m(^rad) with an assumed power-law grain-size distribution n(s) oc s~ 3 ' 5 , which is truncated 
for a minimum and maximum range of s m i n = 0.1 fim and s max = 1mm, respectively (see 
Appendix 1X1). 

For illustrative models, we adopt the standard surface- density distribution of the MMSN 
model. For the gas-to-dust ratio /a, we consider two limiting cases: 1) a constant ratio 
/d = 0.014 and 2) a ratio which includes the effect of ice sublimation. In the latter case, 
f d = 0.0043 for T > 160 K and f d = 0.014 for T < 160 K. In the actual implementation, 
these values are smoothly connected around the sublimation temperature by a hyperbolic 
tangent function. Note that even in the constant-/d cases, the opacity undergoes a transition 
across T = 160 K. 

We perform a set of calculations with this realistic opacity. Similar to the previous 
models with constant opacity, the disk evolves into a quasi-periodic state after about 10 
times of the thermal timescale. First, we show the results for the case with a constant gas- 
to-dust ratio. Figure [12] shows the time evolution of T m in the asymptotic quasi-periodic 
state of the disk with M — 1 x 10 -8 M Q yr _1 . The surface temperature T s is also shown in 
Figure O From equations (jSj) and ffTT]) . we note that T s is independent of time. Comparing 
with the constant-opacity case in which /3 — (Fig. [H]), T s is everywhere larger for the 
realistic grain opacity, because of the effect of superheating. There is a small jump at the 
radius where opacity law changes (T s ~ 160 K). For the same value of M, the structure of 
wave propagation region is very similar to that of the constant-opacity case. 

With a sufficiently large mass-accretion rate (M = 1 x lO _7 M yr -1 ), the disk with a 
realistic opacity is also stabilized by the effect of viscous dissipation (Fig. [13]). There are 
three local maxima in the T m distribution. The innermost peak at about 6 AU corresponds 
to an opacity transition in the disk's surface layer (T s ~ 160 K). All these local peaks do not 
propagate but fluctuate quasi-periodically with a small amplitude. 

Finally, we show the effect of ice sublimation. Figure [141 displays the evolution of T m for 
the model in which the dust surface density is modified by the ice sublimation. Taking into 
account the effect of dust's size sorting on the opacity, we assume that surface dust grains 
has a smaller maximum size s Sjm i n than those in disk interior s m m i n = s m ; n = 1 mm. We 
choose s s ,min = 1 /^ m - Owing to the increase of A s associated with the ice condensation in 
surface layer, the T m distribution attains a local maximum near 5AU. This peak does not 
propagate in time. Interior to this snow line, a wave propagates to about 0.5 AU. Outside 
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this snow line, there is another local maximum near 20 AU, which may be induced by the 
emergence of the first peak. Compare to the model in which the ice sublimation is neglected 
(Fig. [T2l . the ice-condensation induces the formation of a local thermal maximum which 
prevents waves from emerging at large radii and propagate inward. 

5. Summary and Discussion 

We have performed a set of radial one-dimensional calculations to examine the thermal 
evolution of hydrostatic disks, using the direct integration of optical depths t s (T*) to deter- 
mine the optical surface z s and total emitting area-filling factor A s of a superheated layer. 
Our results suggest that, in regions with modest and steep radial temperature gradients, 
the constant x = z s /h assumption is incompatible with the computed height of the surface 
where t s (T*; r, z s ) = 1. The initial state obtained by a fixed- x iteration evolves spontaneously 
to the state where thermal waves grows. The disks evolve to a quasi-periodic state where 
thermal waves continuously propagate toward the star in the intermediate radii. 

The driving mechanism for this thermal instability is the intense stellar irradiation high 
up in the disk atmosphere. It is a consequence of a "shadowing effect" in which the surface 
where most of the stellar photons are intercepted at any given radial location may be affected 
by the vertical structure in the disk regions interior to that radius. This quasi-periodic state 
is stabilized by viscous dissipation associated with the mass accretion flow through the disk. 
For the cases of M — 1O~ 7 M yr _1 wave excitation and propagation are suppressed and the 
disks reach approximately steady states. In order to eliminate the possibility of artificial 
dependence on the initial conditions, we perform the following numerical experiments. By 
setting the initial condition for an approximately steady state with M = 1O~ 7 M yr -1 and 
by decreasing the mass accretion rate to M — 1O~ 8 M yr _1 , we calculate the time evolution 
of the disk. We find that the system evolves to a quasi-periodic state within about 10 times 
of thermal timescale. The asymptotic quasi-periodic state is identical to that obtained with 
the standard calculation for M = lO~ 8 M yr _1 , which is shown in Figure M Inversely, we 
also start a calculation from a quasi-periodic state with M = 1O~ 8 M yr -1 and increase 
the mass accretion rate to M = lO _7 M yr~ 1 . The disk evolves into an approximately 
steady state, which is almost identical to that obtained with the standard calculation for 
M = lO _7 M yr _1 , which is shown in Figure [TUJ These results show that the state of the 
disks are determined by the instantaneous mass accretion rate (note that viscous evolution 
time is much longer than thermal timescale). 

The result that the quasi-periodic wave-propagating states exist in disks with modest 
disk accretion rate is robust to variations in the disks' structural parameters such as their 
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surface density profile, opacity law, and vertical dust distribution. Whereas these parame- 
ters weakly affect the positions of the inner and outer boundaries of the wave-propagating 
domains, they do not affect the basic features of the waves such as their amplitude and 
shape. 

We also perform some calculations for the constant-opacity case using a general surface 
density profile (eq. |B6j instead of eq. [IS]). The results are quite similar with that of the 
simple density profile. The difference is even less noticeable than that due to change of 
opacity law. 

Using the calculated thermal structure of the disk, we obtain the disk SED. We assume 
that the disk is face-on and truncated both at inner 0.1 AU and outer 100 AU. Figure [151 
shows the evolution of the SEDs associated with the model illustrated in Figure dH Contri- 
bution from the central star to the SEDs is included in this figure. The SEDs are expected 
to oscillate periodically within mid-infrared wavelengths. This variation correspond to the 
periodic propagation of the thermal waves. In this case, the period of oscillation is about 
i th0 — 53 yr. In general this time scale varies depending on the surface density distribution 
and accretion rate in the disk as well as the luminosity of the central star. The change of disk 
SED comes from not only emission from the disk interior where T m changes but also emission 
from the superheated surface layer where, even if T s is kept constant, A s would modulate. 
The surface emission also contributes to the water-ice and silicate emission bands. We ex- 
pect the relative heights of these emission bands in SEDs to modulate with time. Note that 
this calculated SED is somewhat artificial because disk is truncated at the inner and outer 
edges. Thus, the direct comparison with observation may be meaningless. Nevertheless, the 
predicted relative variations of SEDs due to propagation of the thermal waves are likely to 
be important in the interpretation of the observed SEDs, especially th ose of the so-calle d 



transitional disks taken by SPITZER Infrared Spectrograph (IRS) (e.g.. iFurlan et al.ll2006l ). 

In this work, we assume t t h to be much longer than ^l^ 1 , so that a hydrostatic equilibrium 
would be quickly re-established after the passage of the thermal waves. This assumption is 
invalid in the outer regions (say > 20 AU, see eq. [7]). But our calculations show that the 
wave propagating region is extended to regions interior to 1 AU. This result implies that 
even if dynamic effect may suppress the thermal waves in the outer d isk, the waves can still 



be excited in the inner regions. According to previous linear analysis (iDullemondl 120001 ) . the 
disk becomes unstable to infinitesimal hydrodynamic perturbations when t t h -C ^k 1 - This 
type of instabilities may also excite the thermal waves. 

The results presented here is based on simple radial one-dimensional analysis in which 
a two-temperature approximation is adopted to describe the vertical structure of the disk. 
We also neglected modulations in the surface density and accretion rate throughout the 
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disk. Our next task is to relax these assumptions and to generalize our results to a set 
of genuine two-dimensional simulations in which the radiation and mass transfer can be 
considered simultaneously. Since dust growth time or radial migration time are comparable 
to the thermal timescale, it will als o be important to consider th e evolution of the dust 
particles (e.g.. iTakeuchi fc Lmll2003l ; iDullemond fc Dominikl l2004bl ). These investigations 
will be carried out in the future and presented elsewhere. 
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A. Opacities 

The Planck mean opacity is given by 

J B v (T TaA )dv 

where subscript j represents "s" (surface) or "m" (disk interior), v is the frequency, B v {T r3i &) 
is the Planck function, and 



f Smax K V (T~, s)s 3 n(s) ds 

MTj) = Js - r l J 3 \,y ■ (A2) 
/ s 6 n{s) ds 

^ ^min 

Here n(s) ds is the number density of grains with radii between s and s + ds, and the size 
distribution has the lower cutoff s min and the upper cutoff s max - In addition, k v (T, s) is the 
single-sized (s) monochromatic (z/) dust opacity of grains with temperature T. 

We assume a power-law size distribution n(s) oc s~ 3,5 , which is truncated for a minimum 
and maximum range of s m i n = 0.1 /xm and s max = 1 mm. 



The grain compositions and optical constants are adopted from iTanaka et al.l (120051 ) 
and references therein. We use the resultant opacity table given by H. Tanaka. The table 
gives two single-sized monochromatic dust opacities as functions of v and s: one for grains 
without ice or organics at high temperatures of T > 160 K and the other for grains including 
ice or organics at T < 160 K. 
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B. Exact Solution of Dust Density Distribution 



The gas density distribution p g of the two-temperature model is given by 

o g (r, 0) exp (-5pJ if M < 



p g (r,z) 



p g (r,z s jexp 



2H 2 



if |Z| > Zg, 



:bv. 



where h = c m fi K (see eq. [5j) and = Cg^ 1 (c s is the sound speed in the surface layer) are 
the gas scale height in disk interior and in the surface layer, respectively. 

In a steady state, the sedimentation flux of the dust grains with their terminal velocity 
balances the diffusive mass flux due to gas turbulence (see eq. [30] in lTakeuchi &: Linll2002l ) 
such that: 

- p d £l K t stop z = — — , (B2) 

where 4to P is the stopping time normalized by f^ 1 , v is the turbulent viscosity, and Sc is the 
Schmidt number representing coupling strength between grains and gas. The nondimensional 
stopping time £ stop is given by 



t 



stop 



PgCt 



where p mat is the material mass density, s is the dust radius, and c t is the mean thermal 
velocity. 



Solving equation (1B2j) with equations (IB II) (in the case of \z\ < z ^) and (|B3D. we obtain 
the dust density distribution of the disk interior to be (see eq. [31] in lTakeuchi fc Linll2002l ) 



Pdm(r, z) = p dm (r, 0)exp 



Set 



stop,m 



2h 2 



a 



exp 



2h 2 



(B4) 



where t s to P ,m is the stopping time in the midplane and a = v/(c m h) is the no n-dimensional 
turbulent viscosity in the so-called a-prescription (jShakura fc Sunyaevlll973l ). The dimen- 
sionless stopping time in the midplane is given by 

7T Pmat s 
2 S 



t 



(B5) 



stop,m 

The midplane dust density Pdm(^, 0) is determined by the vertical integration of equation (IB4j) 
to be the surface density of dust Ed = /dS, where fd is the dust fraction in surface density. 

Solving equation (IB2p with equations (IB1I) (in the case of \z\ > z s ) and (1B3I) . we obtain 
the dust density distribution in the surface layer 



Pds(r, z) = pd m {r,z s )exp 



2H 2 



Set 



stop,m 



h 



aH 



eXP (2^) 1 6XP 



2H 2 



1 



(B6) 
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where we put pds(r, z s ) = pdm(f, z B ) and assume that Sc and a in the surface layer have the 
same values as those in the disk interior. 

In the disk interior, dust sedimentation is not so important unless the radii of dust are 
not so large that we can formally set pdm = Pa- I n the surface layer, dust settling reduces the 
dust density while high surface temperature T s raise it, so that pd also gives a good estimate. 
Thus, instead of equations (1B4I) and (1B6I) . we adopt, in most of the calculations, the simple 
density distribution given by equation (TTST) . 



C. The validity of the grazing-angle approximation 

We derive the grazing-angle approximation (eq. [6]) and discuss its validity. Substituting 
Pd(f, z) in equation (|T5l) with that in equation (fTHI) and assuming A s <C 1, we obtain 

A = T - erfc fe) = T " erfc (^) • (C1) 

where r v = K s (T*)£d/2 and x — z s/h. 

We derive the relation of A s with a path integration from the star to the point (r, z). 
Using equations (fTUl) and (|T8l) . equation ([9]) can be written as 

2 ,i/ 2 /2 r £' c 



r s (T,;r,^(r)) = r v (l + C s 2 ) 1/ %/- / T^e"*^ dr ' = 1 (C2) 



^ JR* K r ') 

with £ s = z 8 /r, S d = Sd(r')/Sd(r), and x' — (s r> ' /h{r'). Here we assume that k s (T*) is 
independent of position of the path. Substituting the integral variable from r' to x' and 
noting that x' ^> 1 for r' = R+, we obtain 

where we define (h = h/r and (' h = h(r')/r'. Here we assume x' is a monotonically increasing 
function with r'. From equations ( 1C1I) and (1C3f) . we obtain 

- ^CiVY' (C4) 



9(X').X 

with 



= [erfc (x/^)] v^exp {- X ' 2 /2) , (C5) 
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where (X) fi X ') tX = X(x')f(x') dx' represents the weighted average with a weight function 



If we substitute the Gaussian weight g(x') with the delta function S(x' — x) an d assume 
Cs <C 1, we obtain 



Except for a small difference, this equation corresponds to equation (Q. However, if the 
gradient dln(h/dlnr changes rapidly in r direction, the approximation used to derive equa- 
tion flC6p is no longer justified. In this case A s is determined not only by the local gradient 
of (h but by the gradients in inner regions because the weight function g(x') is extended to 
the inner radii. For this reason we use equations ([9]) and (|T5|) instead of equation (1C6j) . 
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T m /K 

Fig. I. — Disk thermal emission F m (solid line) and intercepted stellar fluxes F s at r = 10 AU 
with constant surface heights (dashed curves) of z s /r = 0.12, 0.13, and 0.14 (left to right) 
and with constant x — z s/h (dot-dashed line) as functions of disk temperature T m . All fluxes 
are normalized by oT^. 
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Fig. 2. — Initial stage of evolution of the midplane temperature T m for the case of constant 
opacity and (3 — 0. The thick curves cover the range i = i/ith,o — 0.0-1.6 with the interval 
At = 0.2. The surface temperature T s is represented by a thin line in the back panel. The 
surface density distribution is that prescribed by the MMSN model in which p = 1.5. The 
viscous dissipation associated with mass accretion is neglected (i.e., M — 0). 
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Fig. 3. — Time evolution of T m after it has reached a quasi-periodic state. This model is 
the continuation of that shown in Fig. [2] to an epoch i = 8.0-9.6. Each curve is separated 
by At = 0.2. The T s distribution is also shown in a thin line in the back panel. Other 
parameters are same as Fig. [2j 




Fig. 4. — Time evolution of ( s = z s /r in the quasi-periodic state at the same epoch as Fig. [3j 




Fig. 5. — Time evolution of % = z s /h in the quasi-periodic state at the same epoch as Fig. [3j 




Fig. 6. — Time evolution of surface filling factor A s in quasi-periodic state at the same epoch 
as Fig. [31 




Fig. 7. — Time evolution of the logarithmic pressure gradient d In P/d In r in quasi-periodic 
state at the same epoch as Fig. [3j 
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Fig. 9. — Quasi-periodic evolution of T m for steady accretion disk with M = 1O~ 8 M yr _1 
(at i = 14.4-16.0 with At = 0.2). The distribution of T s is also shown with a thin solid line 
in the back panel. Other parameters are same as Fig. [3j 
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Fig. 10. — Quasi-periodic evolution of T m for steady accretion disk with M = 1O _7 M yr _1 
( t = 11.2-14.4 with At = 0.4). The distribution T s is also shown with a thin solid line in 
the back panel. Other parameters are same as Fig. [3j 
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Fig. 11. — Quasi-periodic evolution of T m for a steady accretion disk with p = 1.0 and 
M = lO^Moyr" 1 (t = 9.6-12.8 with At = 0.4). The distribution of T s is also shown with 
a thin solid line in the back panel. Opacities are same as Fig. |3j 
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Fig. 12. — Quasi-periodic evolution of T m for realistic opacities as well as T s (a thin solid line 
in the back panel). The disk structural parameters is that of a MMSN model (p = 1.5) with 
steady M = 10 _8 M Q yr" 1 . The opacity law changes at 160 K but the effect of ice sublimation 
is neglected. Maximum size of surface dust grains is assumed to be 1 mm (same as that in 
disk interior). Times are same as Fig. [3l 




Fig. 13.— Same as Fig. E21 but for a disk with M = 1(T 7 M yr" 1 and i = 8.0-11.2 with 
At = 0.4. 




Fig. 14. — Same as Fig. [121 but taking into account the effect of ice sublimation. The 
maximum size of surface dust grains is reduced to 1 (Ltm. 
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Fig. 15. — Quasi-periodic change of SED calculated from the model shown in Fig. [TU 
Different lines represent the epochs i = 8.4 (thick solid curve), 8.6 (thin solid curve), 8.8 
(dashed curve), 9.0 (dash-dotted curve), 9.2 (dotted curve), and 9.4 (thick solid curve, coincide 
with that at t = 8.4), respectively. Note that t = 53tyr. 



